surge analysis on IOC tide gauge¶
case study: cyclone FUNG WONG
Install¶
Libraries needed for this notebook:
pip install searvey hvplot utide ipykernel geoviews pyarrow
In [31]:
import searvey
import utide
import pandas as pd
import hvplot.pandas
from shapely.geometry import box
In [33]:
ioc_stations = searvey.get_ioc_stations()
ioc_stations
Out[33]:
| ioc_code | gloss_id | lat | lon | country | location | connection | contacts | dcp_id | last_observation_level | ... | observations_expected_per_week | observations_ratio_per_week | observations_arrived_per_month | observations_expected_per_month | observations_ratio_per_month | observations_ratio_per_day | sample_interval | average_delay_per_day | transmit_interval | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | aarh | <NA> | 56.150 | 10.220 | Denmark | Aarhus | ftp | Danish Meteorological Institute ( Denmark ) | NaN | 0.26 | ... | 1008.0 | 99 | 4407 | 4464.0 | 99 | 97 | 10' | 6' | 10' | POINT (10.22 56.15) |
| 1 | abas | 327 | 44.020 | 144.290 | Japan | Abashiri | SWJP40 | Japan Meteorological Agency ( Japan ) | ABASHIRI | 2.10 | ... | 10080.0 | 100 | 43970 | 44640.0 | 98 | 98 | 1' | 8' | 10' | POINT (144.29 44.02) |
| 2 | abed | <NA> | 57.140 | -2.080 | UK | Aberdeen | ftp | National Oceanography Centre ( UK ) | NaN | 3.25 | ... | 672.0 | 100 | 2975 | 2976.0 | 100 | 100 | 15' | 28' | 15' | POINT (-2.08 57.14) |
| 3 | abur | 82 | 31.580 | 131.410 | Japan | Aburatsu | SWJP40 | Japan Meteorological Agency ( Japan ) | ABURATSU | 2.50 | ... | 10080.0 | 100 | 43970 | 44640.0 | 98 | 98 | 1' | 8' | 10' | POINT (131.41 31.58) |
| 4 | acaj | 182 | 13.574 | -89.838 | El Salvador | Acajutla SV | SZXX01 | Ministerio de Medio Ambiente y Recursos Natura... | 300434064008810 | 1.16 | ... | 10080.0 | 91 | 36760 | 44640.0 | 82 | 90 | 1' | 6' | 5' | POINT (-89.838 13.574) |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1720 | zkth | <NA> | 37.781 | 20.905 | Greece | Zakynthos, Ionian Islands | bgan | National Observatory of Athens ( Greece ) | GR-ZKTH-00 | 0.01 | ... | 10080.0 | 0 | NaN | NaN | 0 | 0 | 1' | NaN | 1' | POINT (20.905 37.781) |
| 1721 | zldw | <NA> | 51.346 | 3.200 | Belgium | Zeebrugge Leopold II dam | web | Maritieme Dienstverlening en Kust ( Belgium ) | NaN | 3.86 | ... | 2016.0 | 100 | 8902 | 8928.0 | 100 | 100 | 5' | 8' | 5' | POINT (3.2 51.346) |
| 1722 | zwdw | <NA> | 51.355 | 3.178 | Belgium | Zeebrugge Wielingendok | web | Maritieme Dienstverlening en Kust ( Belgium ) | NaN | 3.88 | ... | 2016.0 | 100 | 8908 | 8928.0 | 100 | 100 | 5' | 8' | 5' | POINT (3.178 51.355) |
| 1723 | zygi | <NA> | 34.727 | 33.338 | Cyprus | Zygi | ftp | Cyprus Oceanography Center ( Cyprus ) | NaN | 1.97 | ... | 20160.0 | 0 | -down- | NaN | 0 | 0 | 0.5' | NaN | 1' | POINT (33.338 34.727) |
| 1724 | zygi1 | <NA> | 34.726 | 33.340 | Cyprus | Zygi | bgan | Cyprus Marine and Maritime Institute ( Cyprus ... | ZYGI1 | 0.63 | ... | 10080.0 | 0 | -down- | NaN | 0 | 0 | 1' | NaN | 1' | POINT (33.34 34.726) |
1725 rows × 25 columns
In [34]:
bbox = box(100.0, -5.0, 130.0, 25.0) # lon_min, lat_min, lon_max, lat_max
In [36]:
selected_stations = ioc_stations[ioc_stations.within(bbox)]
selected_stations
Out[36]:
| ioc_code | gloss_id | lat | lon | country | location | connection | contacts | dcp_id | last_observation_level | ... | observations_expected_per_week | observations_ratio_per_week | observations_arrived_per_month | observations_expected_per_month | observations_ratio_per_month | observations_ratio_per_day | sample_interval | average_delay_per_day | transmit_interval | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 45 | ambon | 68 | -3.683 | 128.183 | Indonesia | Ambon | SZXX01 | Geospacial Agency of Indonesia ( Indonesia ) +... | 301434060409970 | 4.46 | ... | 10080.0 | 91 | 36510 | 44640.0 | 82 | 90 | 1' | 5' | 5' | POINT (128.183 -3.683) |
| 47 | AMTSI | <NA> | -0.847 | 121.601 | Indonesia | Port of Ampana, Central Sulawesi | web | Kepala Badan Meteorologi, Klimatologi, dan Geo... | NaN | 0.01 | ... | 10080.0 | 96 | 37857 | 44640.0 | 85 | 98 | 1' | 2' | 5' | POINT (121.601 -0.847) |
| 132 | BETSI | <NA> | 4.230 | 126.788 | Indonesia | Port of Beo, North Sulawesi | web | Kepala Badan Meteorologi, Klimatologi, dan Geo... | NaN | 0.01 | ... | 10080.0 | 96 | 43470 | 44640.0 | 97 | 76 | 1' | 1' | 5' | POINT (126.788 4.23) |
| 136 | bitu | 69 | 1.439 | 125.190 | Indonesia | Bitung | SZXX01 | Geospacial Agency of Indonesia ( Indonesia ) +... | 300434067058990 | 8.15 | ... | 10080.0 | 90 | 36675 | 44640.0 | 82 | 89 | 1' | 5' | 5' | POINT (125.19 1.439) |
| 183 | bupo | <NA> | -1.030 | 100.396 | Indonesia | Bungus Port | bgan | Joint Research Centre ( Europe ) +Ministry of ... | BUPO | 1.45 | ... | 10080.0 | 0 | -down- | NaN | 0 | 0 | 1' | NaN | 1' | POINT (100.396 -1.03) |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1621 | txil | <NA> | 22.353 | 120.383 | Taiwan | Xiao Liuqiu | web | Central Weather Administration ( Taiwan ) | NaN | 0.20 | ... | 168.0 | 100 | 743 | 744.0 | 100 | 100 | 60' | 16' | 1h | POINT (120.383 22.353) |
| 1622 | tyon | <NA> | 22.819 | 120.198 | Taiwan | Yongan | web | Central Weather Administration ( Taiwan ) | NaN | -999.00 | ... | 168.0 | 100 | 743 | 744.0 | 100 | 100 | 60' | 16' | 1h | POINT (120.198 22.819) |
| 1672 | vung | <NA> | 10.340 | 107.071 | Viet Nam | Vung Tau | SZXX01 | Hydro-meteorological and Environmental Station... | 300434067059990 | 2.26 | ... | 10080.0 | 88 | 35835 | 44640.0 | 80 | 87 | 1' | 5' | 6' | POINT (107.071 10.34) |
| 1706 | xish | <NA> | 16.830 | 112.330 | China | Xisha | SZCI01 | China Meteorological Administration ( China ) | 11745 | 1.30 | ... | 10080.0 | 0 | -down- | NaN | 0 | 0 | 1' | NaN | 5' | POINT (112.33 16.83) |
| 1717 | zhap | 78 | 21.580 | 111.820 | China | Zhapo | SZCI01 | China Meteorological Administration ( China ) | 09731 | 2.42 | ... | 10080.0 | 44 | 18491 | 44640.0 | 41 | 45 | 1' | 7' | 1' | POINT (111.82 21.58) |
70 rows × 25 columns
In [32]:
plot_ = selected_stations.hvplot.points(
geo=True,
tiles=True,
hover_cols=["ioc_code"],
s=100,
line_color='k',
title = "stations in the selected region"
).opts(width=600, height=700)
plot_
Out[32]:
detide the stations¶
In [30]:
def surge(ts: pd.Series, lat: float, rsmp: int = None):
ts0 = ts.copy()
OPTS = {
"constit": "auto",
"method": "ols", # ols is faster and good for missing data (Ponchaut et al., 2001)
"order_constit": "frequency",
"Rayleigh_min": 0.97,
"lat": lat,
"verbose": True,
}
if rsmp is not None:
ts = ts.resample(f"{rsmp}min").mean()
ts = ts.shift(freq=f"{rsmp / 2}min")
coef = utide.solve(ts.index, ts, **OPTS)
tidal = utide.reconstruct(ts0.index, coef, verbose = OPTS['verbose'])
return pd.Series(data=ts0.values - tidal.h, index=ts0.index)
download the data¶
create a data folder with raw and surge subfolders
In [12]:
! mkdir -p data/{raw,surge}
let's first download all the station in the bounding box
In [16]:
drop_columns = ["bat"]
def serialize(d): # to export metadata
out = {}
for k, v in d.items():
if v is not pd.NA and not pd.isna(v):
out[k] = str(v)
return out
In [17]:
for irow, row in selected_stations.iterrows():
lat = row.lat
ioc_code = row.ioc_code
if ioc_code in list(detide_selection.keys()):
df_raw = searvey.fetch_ioc_station(
ioc_code,
pd.Timestamp.now()-pd.Timedelta(days=90), # we need at least 90 days (in theory we'd need more..) to remove properly tidal constituents
pd.Timestamp.now()
)
df_raw.attrs = {**serialize(dict(row)), **{"signal_type": "raw"}}
df_raw.to_parquet(f"data/raw/{ioc_code}.parquet")
In [6]:
for irow, row in selected_stations.iterrows():
df_raw = pd.read_parquet(f"data/raw/{row.ioc_code}.parquet")
df_raw.loc["2025-11-01":].dropna().hvplot(
title = f"detided signal for {row.ioc_code}, {row.location}, lat:{row.lat}, lon:{row.lon}"
)
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay [Variable]
Out[6]:
detide the station¶
Let's look at the stations of interest
In [ ]:
In [18]:
detide_selection = {
"lega" :"rad",
"luba" :"prs",
"quin" :"ras",
"qing" :"rad",
"quar" :"flt",
"mani" :"prs",
"subi" :"rad",
"thsi" :"rad",
"txil" :"rad",}
In [27]:
for irow, row in selected_stations.iterrows():
if row.ioc_code in list(detide_selection.keys()):
df_raw = pd.read_parquet(f"data/raw/{row.ioc_code}.parquet")
df_surge = surge(df_raw[detide_selection[row.ioc_code]].dropna(), lat=row.lat, rsmp=2)
df_surge.loc["2025-11-01":].hvplot(
title = f"detided signal for {row.ioc_code}, {row.location}, lat:{row.lat}, lon:{row.lon}"
)
df_surge.attrs = dict(row)
df_surge.attrs = {**serialize(dict(row)), **{"signal_type": "detide"}}
df_surge.to_frame().to_parquet(f"data/surge/{row.ioc_code}.parquet")
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done. prep/calcs ... done.
Out[27]:
In [26]:
(plot_*ioc_stations[ioc_stations.ioc_code.isin(list(detide_selection.keys()))].hvplot.points(
geo=True,
tiles=True,
c = 'r',
s=50,
hover_cols=["ioc_code"],
title = "stations that recorded a surge"
)).opts(width=800, height=800)
Out[26]:
look at a specific station¶
In [ ]:
station = "lega"
for irow, row in selected_stations.iterrows():
if row.ioc_code == station:
df_raw = pd.read_parquet(f"data/raw/{row.ioc_code}.parquet")
df_raw
station
df_raw.loc["2025--01":].dropna().hvplot(
title = f"signal for {row.ioc_code}, {row.location}, lat:{row.lat}, lon:{row.lon}"
)